Life and Death
This simulation actually extends an example from the DifferentialEquations.jl documentation describing a growing cell population. Here we take those underlying model dynamics and combine them with CellularPotts.jl
using CellularPotts, DifferentialEquationsCellularPotts.jl Setup
Start by defining a space and the characteristics of the cells
space = CellSpace(200,200)
initialCellState = CellTable(
[:Epithelial],
[200],
[1]);The single cell will be positioned at the halfway point within the space.
positions = [size(space) .÷ 2]
initialCellState = addcellproperty(initialCellState, :positions, positions)┌────────────┬─────────┬─────────┬─────────┬────────────────┬────────────┬───────────────────┬─────────────────────┐
│ names │ cellIDs │ typeIDs │ volumes │ desiredVolumes │ perimeters │ desiredPerimeters │ positions │
│ Symbol │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Tuple{Int64, Int64} │
├────────────┼─────────┼─────────┼─────────┼────────────────┼────────────┼───────────────────┼─────────────────────┤
│ Medium │ 0 │ 0 │ 0 │ 0 │ 0 │ 0 │ (100, 100) │
│ Epithelial │ 1 │ 1 │ 0 │ 200 │ 0 │ 168 │ (100, 100) │
└────────────┴─────────┴─────────┴─────────┴────────────────┴────────────┴───────────────────┴─────────────────────┘
As per the growing cell population example, we define a theoretical protein X that increases linearly in time with rate parameter α
const α = 0.3;Include an initial protein concentration for :Medium and set it to zero
u0 = [0.0, 0.2]
initialCellState = addcellproperty(initialCellState, :ProteinX, u0)┌────────────┬─────────┬─────────┬─────────┬────────────────┬────────────┬───────────────────┬─────────────────────┬──────────┐
│ names │ cellIDs │ typeIDs │ volumes │ desiredVolumes │ perimeters │ desiredPerimeters │ positions │ ProteinX │
│ Symbol │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Tuple{Int64, Int64} │ Float64 │
├────────────┼─────────┼─────────┼─────────┼────────────────┼────────────┼───────────────────┼─────────────────────┼──────────┤
│ Medium │ 0 │ 0 │ 0 │ 0 │ 0 │ 0 │ (100, 100) │ 0.0 │
│ Epithelial │ 1 │ 1 │ 0 │ 200 │ 0 │ 168 │ (100, 100) │ 0.2 │
└────────────┴─────────┴─────────┴─────────┴────────────────┴────────────┴───────────────────┴─────────────────────┴──────────┘
Keeping the model simple, we'll only include an Adhesion and Volume penalty
penalties = [
AdhesionPenalty([0 30;
30 30]),
VolumePenalty([5]),
]2-element Vector{Penalty}:
AdhesionPenalty([0 30; 30 30])
VolumePenalty([0, 5])Declare the model and position the cells in the model
cpm = CellPotts(space, initialCellState, penalties);
positionCells!(cpm)DifferentialEquations.jl setup
Use a periodic callback to increment the CPM model every time step
Currently there isn't a simple method to log states in CellPotts models (work in progress! 🙂). For now, we need to create an extrernal variable to log how the nodeIDs change over time.
spaceLog = [cpm.space.nodeIDs];Next we need to step the CPM in time at regular intervals. We can accomplish this during the differential equations solve using a periodic callback. This function will trigger at regular interval according to the time scale. Larger timescales correspond to faster cell movement.
Step the model and update the space log
function cpmUpdate!(integrator, cpm, spaceLog)
ModelStep!(cpm)
push!(spaceLog, copy(cpm.space.nodeIDs))
endcpmUpdate! (generic function with 1 method)Set the timescale and create the callback
timeScale = 100
pcb = PeriodicCallback(integrator -> cpmUpdate!(integrator, cpm, spaceLog), 1/timeScale)SciMLBase.DiscreteCallback{DiffEqCallbacks.var"#44#48"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.PeriodicCallbackAffect{Main.##4812.var"#1#2", Float64, Base.RefValue{Float64}, Base.RefValue{Int64}}, DiffEqCallbacks.var"#45#49"{Bool, DiffEqCallbacks.var"#46#50"{Bool}, Float64, DiffEqCallbacks.PeriodicCallbackAffect{Main.##4812.var"#1#2", Float64, Base.RefValue{Float64}, Base.RefValue{Int64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(SciMLBase.FINALIZE_DEFAULT)}(DiffEqCallbacks.var"#44#48"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}(0.01, Base.RefValue{Int64}(0), Base.RefValue{Float64}(Inf)), DiffEqCallbacks.PeriodicCallbackAffect{Main.##4812.var"#1#2", Float64, Base.RefValue{Float64}, Base.RefValue{Int64}}(Main.##4812.var"#1#2"(), 0.01, Base.RefValue{Float64}(Inf), Base.RefValue{Int64}(0)), DiffEqCallbacks.var"#45#49"{Bool, DiffEqCallbacks.var"#46#50"{Bool}, Float64, DiffEqCallbacks.PeriodicCallbackAffect{Main.##4812.var"#1#2", Float64, Base.RefValue{Float64}, Base.RefValue{Int64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}(false, DiffEqCallbacks.var"#46#50"{Bool}(false), 0.01, DiffEqCallbacks.PeriodicCallbackAffect{Main.##4812.var"#1#2", Float64, Base.RefValue{Float64}, Base.RefValue{Int64}}(Main.##4812.var"#1#2"(), 0.01, Base.RefValue{Float64}(Inf), Base.RefValue{Int64}(0)), Base.RefValue{Int64}(0), Base.RefValue{Float64}(Inf)), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1])This is taken directly from the example. Each cell is given the following differential equation
\[ \frac{\mathrm{d} X}{\mathrm{d} t} = \alpha X\]
function f(du,u,p,t)
for i in eachindex(u)
du[i] = α*u[i]
end
endf (generic function with 1 method)Also from the differential equations example. This callback is tiggered whenever a cell's ProteinX is greater than 1
condition(u,t,integrator) = 1-maximum(u)
function affect!(integrator,cpm)
u = integrator.u
resize!(integrator,length(u)+1)
cellID = findmax(u)[2]
Θ = rand()
u[cellID] = Θ
u[end] = 1-Θ
#Divide the cells in the CPM
CellDivision!(cpm, cellID-1)
return nothing
end
ccb = ContinuousCallback(condition,integrator -> affect!(integrator, cpm));Collect the callbacks together
callbacks = CallbackSet(pcb, ccb);Define the ODE model and solve
tspan = (0.0,20.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(),callback=callbacks);We can replicate the plots from the original example
using Plots, PrintfPlot the total cell count over time
plot(sol.t,map((x)->length(x),sol[:]),lw=3,
ylabel="Number of Cells",xlabel="Time")
Plot ProteinX dynamics for a specific cell
ts = range(0, stop=20, length=100)
plot(ts,map((x)->x[2],sol.(ts)),lw=3, ylabel="Amount of X in Cell 1",xlabel="Time")
Finally, we can create an animation of the CPM to see the cells dividing. I've dropped the first few frames because the first cell takes a while to divide.
anim = @animate for t in Iterators.drop(eachindex(spaceLog),5*timeScale)
currTime = @sprintf "Time: %.2f" t/timeScale
heatmap(spaceLog[t], axis=nothing,legend = :none, size=(500,500),title=currTime)
end
gif(anim, "LifeAndDeath.gif", fps = 30)This page was generated using Literate.jl.